clear;
close all;

dynare model_file;
NumberOfEndogenousVariables = M_.endo_nbr;                  
for ii = 1:NumberOfEndogenousVariables                        
  varname =   deblank(M_.endo_names(ii,:));               
  temp = strcat('o.', varname, ' = ii;  ');          
  eval(temp{1,1});
end

YY = oo_.endo_simul'; % T*N, simulated data of all variables
YY_m = mean(YY); % mean of each variable
YY_sd = std(YY); % standard deviation of each variable


fit = fitlm( [pi(1:end-1)*400,pie(1:end-1)*400],pi(2:end)*400); % core inflation
shock_pi = fit.Residuals.Raw; % core inflation shock

fit = fitlm( [pi(1:end-1)*400, pie(1:end-1)*400], pihnum(2:end)*400);
shock_pih = fit.Residuals.Raw; % headline inflation shock

fit = fitlm( [pi(1:end-1)*400, pie(1:end-1)*400], pie(2:end)*400);
shock_pie = fit.Residuals.Raw; % energy inflation shock

% nominal return
m(m<0) = 0.1; % set negative 
mn = log(m) - pi;  % nominal SDF
mn_corr = 1-(0.07/2)^2/(2*var(mn)); % vol(Delta s) = vol(m*-m)


% Assuming m*,m have the same variance, rho = 1-vol(Delta s)/(2*vol(m))
mn_f = mn_corr*mn + sqrt(var(mn)*(1-mn_corr^2))*randn(options_.periods,1); % nominal foreign SDF (without specifying foreign inflation)
ds = mn_f - mn; 
ds = ds - mean(ds);
rxfx =  exp(ds(2:end))*(mean(exp(yn1)))-exp(yn1(1:end-1)); 
rxdollar = rxfx.*(mean(exp(yn1))>exp(yn1(1:end-1))) - rxfx.*(mean(exp(yn1))<exp(yn1(1:end-1))); % dollar carry

rmn =  log(exp(rm).*exp(pihnum))  ; % nominal excess return
rmnagg =  log(exp(rmagg).*exp(pihnum))  ; % nominal excess return


RR = [ [exp(rmnagg(2:end)) exp(-19*yn19(2:end))./exp(-20*yn20(1:end-1))] - exp(yn1(1:end-1)) ... 
 rxfx  ...
 pe(2:end).*exp(pi(2:end))./fen(1:end-1)-1  ]; % stock, bond (5y), currency, commodity, excess return

inf_cor = zeros(3,1);
infl = [pi pie pihnum];
for i = 1:3
    inf_cor(i,1) = corr(infl(1:end-1,i),infl(2:end,i));
end
table_11a = [oo_.mean(o.dy)*400;YY_sd([o.dy o.dcr o.di o.dce])'*200;oo_.mean([o.rf o.yn1])*400];
table_11b = [oo_.mean([o.pi o.pie o.pihnum])'*400;YY_sd([o.pi o.pie o.pihnum])*400;inf_cor'];
table_11c = [mean(RR,1)'*400 std(RR,1)'*200];
%% Asset pricing tests
test = RR*400;
factor1 = [shock_pi shock_pie];
[beta1,lambda1,se_beta1,t_beta1,lambda_se1,a1,r21,r2_2step1]...
    =cross_section(test,factor1,1);
factor2 = shock_pih;
[beta2,lambda2,se_beta2,t_beta2,lambda_se2,a2,r22,r2_2step2]...
    =cross_section(test,factor2,1);
table_11d = [beta1' beta2'];
table_11e = [lambda1' lambda2];
table_12 = [YY_sd([o.pi o.pie])'*400;mean(RR(:,1))*400;mean(RR(:,2))*400;beta1(1,1:2)';beta1(2,1:2)';lambda1];
%% Stock-bond correlation
RR2 = [exp(rmnagg(2:end)) exp(-19*yn19(2:end))./exp(-20*yn20(1:end-1))]-exp(yn1(1:end-1)); % stock and bond excess return
RR2 = RR2*400; % annualize
table_13 = [corr(RR2(:,1),RR2(:,2));beta1(1,1:2)';beta1(2,1:2)';YY_sd([o.pi o.pie])'*400];


